Based on the replication code for this paper, this workbook aims to explore Reviewer 2s objection:

“All legislators that remain in the dataset between t and t+1 increase in experience 1 unit. Thus, there is no identifying variation to separate the effect of experience from the effect of time shocks in the within-individual specification.

Thus, the coefficients on the dummy variables indicating different levels of legislative experience are not identified in the OLS regression with time- and legislator fixed effects.”

Replication data

From code/replication.R, which starts with the main data file (dcounts_min.Rdata) and makes transformations for three sets of models:

  1. counts per member-year (dcounts_tenure.Rdata/AgencyComm.dta)
  2. ratio per member-year (dcounts_ratio.Rdata/ProportionContact.dta)
  3. counts per district-year (dcounts_per_district.Rdata, DistrictLevel.dta)

The same three datasets are used in Replication.do

load(here::here("data", "dcounts_tenure.Rdata"))
load(here::here("data", "dcounts_ratio.Rdata"))
load(here::here("data", "dcounts_per_district.Rdata"))
load(here::here("data" , "members.Rdata"))


# var names for tables
dcounts_tenure %<>% 
  mutate(Legislator = icpsr,
         Legislator_x_Agency = icpsr_agency,
         Year_x_Agency = agency_year)

dcounts_tenure %>% head() %>% kablebox()
agency icpsr chamber.x year perYear perYear_con perYear_policy congress party party_code state_abbrev district_code nominate.dim1 presidents_party female chair ranking_minority party_leader party_whip majority prestige prestige_chair yearelected state first_cong first_year tenure first second third fourth fifth sixth chamber.y oversight_committee oversight_committee_chair icpsr_agency agency_year max_year survive chamber n Legislator Legislator_x_Agency Year_x_Agency
ABMC 10713 House 2016 0 0 0 114
100 MI 13 -0.658 1 0 0 1 0 0 0 0 0 2014 michigan 89 1965 51 0 0 0 0 0 0 NA NA NA ABMC_10713 ABMC_2016 53 1 House 1 10713 ABMC_10713 ABMC_2016
ABMC 10713 House 2017 0 0 0 115
100 MI 13 -0.658 0 0 0 1 0 0 0 0 0 2016 michigan 89 1965 52 0 0 0 0 0 0 NA NA NA ABMC_10713 ABMC_2017 53 1 House 1 10713 ABMC_10713 ABMC_2017
ABMC 10713 House 2018 0 0 0 115
100 MI 13 -0.658 0 0 0 1 0 0 0 0 0 2016 michigan 89 1965 53 0 0 0 0 0 0 NA NA NA ABMC_10713 ABMC_2018 53 1 House 1 10713 ABMC_10713 ABMC_2018
ABMC 13035 House 2016 0 0 0 114
100 NY 13 -0.514 1 0 0 0 0 0 0 1 0 2014 new york 92 1971 45 0 0 0 0 0 0 NA NA NA ABMC_13035 ABMC_2016 45 1 House 1 13035 ABMC_13035 ABMC_2016
ABMC 14009 Senate 2016 0 0 0 114
200 MS 0 0.287 0 NA 1 0 0 0 1 1 1 2014 mississippi 93 1973 43 0 0 0 0 0 0 NA NA NA ABMC_14009 ABMC_2016 45 1 Senate 1 14009 ABMC_14009 ABMC_2016
ABMC 14009 Senate 2017 0 0 0 115
200 MS 0 0.287 1 NA 1 0 0 0 1 1 1 2016 mississippi 93 1973 44 0 0 0 0 0 0 NA NA NA ABMC_14009 ABMC_2017 45 1 Senate 1 14009 ABMC_14009 ABMC_2017

Reviewer 2s objection is about the models using (1) counts per member-year, so that is the DV I focus on here.

Values to predict at

The most frequent legislator agency pair: McCain x VA with 2159 letters.

The average is 73 letters per legislator per year, a little over 1.29 letters per agency per year.

For predictions, we need an average member-year-agency triad. The first table below shows the members who averaged 73, entered Congress between 2007 and 2011, were observed in our data for at least seven years, and most often submitted exactly one letter to an agency (the closest integer 1.29). For realism, this person would ideally, at some point in our data, become a chair. Of the modal legislators, Phil Roe meets this last criterion. The IRS is a fairly modal agency and one of the agencies to which Phil Roe submitted one letter in 2015, so predictions are at Phil Roe X IRS X 2015. Other legislators, years, and agencies are simply intercept shifts.

modal <- dcounts_tenure %>% 
  group_by(icpsr) %>% 
  filter(first_year > 2006, first_year < 2013, max_year > 6) %>% 
  mutate(Chair = sum(chair, na.rm = T) >1) %>% 
  group_by(icpsr, year) %>%
  mutate(perLeg = sum(perYear)) %>% 
  filter(perYear == 1, perLeg == 73 ) %>% 
  select(perAgency = perYear, perLeg, icpsr, icpsr_agency, agency_year, first_year, max_year, Chair) %>% 
  distinct() %>% 
  ungroup() %>% left_join(members %>% select(bioname, icpsr))



modal %<>% filter(Chair)


# by chair and presidents' party
values <- tidyr::expand(dcounts_tenure,
                        agency = "Treasury_IRS",
                        chair = chair, 
                        ranking_minority    = FALSE, 
                        prestige = TRUE, 
                        first = FALSE,
                        second = FALSE,
                        third = FALSE,
                        fourth = FALSE,
                        fifth = FALSE, 
                        sixth = FALSE,
                        Legislator_x_Agency = "Treasury_IRS_20947",
                        Year_x_Agency = "Treasury_IRS_2015",
                        Legislator = "20947",
                        Year = "2015",
                        majority = TRUE,
                        presidents_party = presidents_party
                   ) %>% 
  drop_na(majority, chair)

# by year 
values_tenure <- tidyr::expand(dcounts_tenure,
                        agency = "Treasury_IRS",
                        chair = chair, 
                        ranking_minority    = FALSE, 
                        prestige = TRUE, 
                        first = first,
                        second = second,
                        third = third,
                        fourth = fourth,
                        fifth = fifth, 
                        sixth = sixth,
                        Legislator_x_Agency = "Treasury_IRS_20947",
                        Year_x_Agency = "Treasury_IRS_2015",
                        Legislator = "20947",
                        Year = "2015",
                        majority = TRUE,
                        presidents_party = TRUE
                   ) %>% 
  drop_na(chair)

Simulated data

Counts per legislator per agency per year

I replace the dependent variable, the count of letters that each legislator writes to each agency each year with a normal distribution with a mean 1 (the true mean in the data is 1.29).

I also introduce a correlation with the calendar year by adding +1 for each year beyond 2007 (the first year in the data.) In this simulation 2017 would have a mean of 11.

Later, I will introduce a correlation with tenure. (See “100 simulations section” below)

n <- nrow(dcounts_tenure)

sim <- dcounts_tenure %>% 
  # REPLACE DV WITH NORMAL DIST
  mutate(perYear = rnorm(n, year-2007)
         # OTHER VARS LEFT ALONE 
         # chair,
         # first,
         # second,
         # third,
         # fourth,
         # fifth,
         # sixth,
         # Year, 
         # Agency,
         # Legislator
  ) 

dcounts_tenure <- sim

Models

Our primary models are a series of difference-in-differences regressions that examine changes that are within legislator and agency pairs.

Estimates from Model 2 (Column 2 below) provides the estimated effects from the difference-in-differences specification.

$Y_{ijt} = {}^{’} {it} + {s = 1}^{6} {s} ({it} = s) + {ij} + {jt} + m_{it} + p_{it} + _{ijt} $

Where \(Y_{ijt}\) represents the number of requests legislator \(i\) makes to agency \(j\) in year \(t\). Our analysis in this section is at the legislator-agency-year level. \(\gamma_{ij}\) is a fixed effect for the legislator-agency pair. This fixed effect accounts for legislators’ characteristics, such as legislators who are more skillful at filling constituency service requests than other legislators. Critically for our research design, this fixed effect also enables us to account for time-invariant constituent demands, ensuring differences in constituent demand do not drive our results. It also accounts for state and districts characteristics, including population, demographics, and local industries that might be particularly likely to request help with specific agencies. This difference-in-difference design ensures that coefficients \({\beta}\) capture variation related to changes in institutional power or experience, not other factors that may vary across districts, legislators, or agencies. The model also accounts for the different periods for which data were available from each agency. \(\delta_{jt}\) is an agency-year fixed effect. This takes into account agency-level shocks that may affect legislator requests.

Assuming that legislators’ trends in the level of requests follow parallel paths, \(\boldsymbol{\beta}\) represents the average effect of changing institutional power on a legislator’s provision of constituency service.

In this same regression we also include indicators for legislators’ first six years in Congress, $ {s = 1}^{6} {s} _{it}$. The effects of interest \(\eta_{1}, \eta_{2} : \eta_{6}\) describe how a legislator’s provision of constituency service at levels of seniority between one and six years differ from legislators who serve beyond six years. We focus on constituency service levels in each of the first six years of a legislator’s tenure to assess how constituency service changes over their initial years in Congress. This design allows us to assess the extent to which new legislators face start-up costs.

Estimation

Fixed effects general linear model fixest::fixest

feglm (perYear ~ 
 chair + ranking_minority + prestige +
 first + second + third + fourth + fifth + sixth + 
 majority + presidents_party | Year_x_Agency + Legislator_x_Agency, 
 cluster = "Legislator",
 data = dcounts_tenure)

Where

perYear = \(Y_{ijt}\) = Letters per legislator \(i\) per agency \(j\) per year \(t\)

chair= \(\textbf{Committee Position}_{it}\) = Legislator \(i\) is a committee chair in year \(t\)

first + second + third + fourth + fifth + sixth = \(\text{I}\left(\text{tenure}_{it} = s\right)\) = Legislator i is in their \(s\) (first, second, third, fourth, fifth, sixth) year of service in calendar year \(t\)

\(\gamma_{ij}\) = a fixed effect for the legislator-agency pair = Legislator_x_Agency

\(\delta_{jt}\) = an agency-year fixed effect = Year_x_Agency


Counts per member per year

Clustering standard errors at the legislator level:

“Clustering on the panel variable produces an estimator of the VCE that is robust to cross-sectional heteroskedasticity and within-panel (serial) correlation that is asymptotically equivalent to that proposed by Arellano (1987)”

Replicated using fixest::feglm here and xtreg in Replication.do. See more about robust (clustered) standard errors with fixed effects are calculated by fixest here and by xtreg here.

# Coef Map
cm = c("chair" = "Committee Chair",
       "ranking_minority" = "Ranking Member",
       "prestige" = "Prestige Committee",
       "first" = "First Year",
       "second" = "Second Year",
       "third" = "Third Year",
       "fourth" = "Fourth Year",
       "fifth" = "Fifth Year",
       "sixth" = "Sixth Year",
       "majority" = "Majority",
       "presidents_party" = "President's Party",
       "Legislator" = "Legislator", 
       "Agency" = "Agency",
       "FE: icpsr_agency" = "Legislator × Agency Fixed Effects",
       "icpsr_agency" = "Legislator × Agency Fixed Effects")

coef_omit = "(Intercept)|majority|presidents_party"
coef_omit = "none"

setFixest_dict(cm)

gof_omit = "R.*|AIC|BIC|Log.*|Std.*"

# table formatting to match stata (sort of)
format_table <- . %>% 
  row_spec(row = 1, bold = T, hline_after = TRUE) %>% 
      kable_styling(font_size = 11) %>% 
                    #full_width = TRUE, 
                    #latex_options = c("repeat_header")) %>%
  str_replace("Num.Obs.", "Observations") %>% 
  str_replace("Std.Errors", "Robust/Clustered Std. Errors") %>% 
  str_replace("FE: Legislator_x_Agency", "Legislator x Agency FE") %>%
  str_replace("FE: Year_x_Agency", "Year x Agency FE") %>% 
  str_replace("FE: Year", "Year Fixed Effects") %>% 
  str_replace("FE: Legislator", "Legislator Fixed Effects") %>%
    str_replace_all("X|✓", "\\\\\\checkmark") %>% 
  # a random midrule appeared in the wrong place 
  #str_remove("\\midrule") %>% 
  #  extract just the table, no caption etc
  str_extract("\\\\begin\\{tabular\\}(\n|.)*tabular\\}") 

Total Letters

Coeficient Plots

  • Figures: figs/m-total-[1:4].png
  1. “cross-sectional”
  2. “diff-in-diff”
  3. “cross-sectional re-elected”
  4. “diff-in-diff re-elected”
testing = F

# paper table 
# 1 
# cross-sectional 
m_total_cross <-feglm (perYear ~ 
                         chair + ranking_minority + prestige + 
                    first + second + third + fourth + fifth + sixth + 
                    majority + presidents_party | Year_x_Agency,
                    cluster = "Legislator", 
           data = dcounts_tenure)

if(testing){modelsummary(m_total_cross, gof_omit = "R.*", coef_omit = coef_omit)}

coefplot(m_total_cross, horiz = T, drop = "(Intercept)") 



# 2
m_total_dnd <-feglm (perYear ~ 
                       chair + ranking_minority + prestige + 
                    first + second + third + fourth + fifth + sixth + 
                    majority + presidents_party | Year_x_Agency + Legislator_x_Agency, 
                    cluster = "Legislator",
                    # ALT vcov = hetero ~ ssc(cluster.adj = TRUE),
           data = dcounts_tenure)

if(testing){modelsummary(m_total_dnd, gof_omit = "R.*", coef_omit = coef_omit)}

coefplot(m_total_dnd, horiz = T) 


# 3
m_total_2nd <-feglm (perYear ~ 
                       chair + ranking_minority + prestige + 
                    first + second + third + fourth + fifth + sixth + 
                    majority + presidents_party | Year_x_Agency + Legislator_x_Agency, 
                    cluster = "Legislator", 
           data = dcounts_tenure %>% filter(survive == 1))

if(testing){modelsummary(m_total_2nd, gof_omit = "R.*", coef_omit = coef_omit)}
coefplot(m_total_2nd, horiz = T) 


# 4
m_logtotal_dnd <-feglm (log(perYear + 1) ~ 
                          chair + ranking_minority + prestige + 
                    first + second + third + fourth + fifth + sixth + 
                    majority + presidents_party | Year_x_Agency + Legislator_x_Agency, 
                    cluster = "Legislator", 
           data = dcounts_tenure)

if(testing){modelsummary(m_logtotal_dnd, gof_omit = "R.*", coef_omit = coef_omit)}

coefplot(m_logtotal_dnd, horiz = T) 
Total Letters Per YearTotal Letters Per YearTotal Letters Per YearTotal Letters Per Year

Total Letters Per Year

Total Letters Table

  • Table: tables/models_total.tex (sans stars)
models_total <- list(
  "(1)" = m_total_cross,
  "(2)" = m_total_dnd,
  "(3)" =  m_total_2nd,
  "(4)" = m_logtotal_dnd
)

rows <- tibble(
  term = c("Dependent Variable", 
           "Majority",
           "President's Party",
           "All Legislators", 
           "Served At Least 2nd Term"),
  `(1)` = c("Count", "✓", "✓", "✓", ""),
  `(2)` =c("Count", "✓", "✓", "✓", ""),
  `(3)` = c("Count", "✓","✓","", "✓"),
  `(4)` = c("Log(Count+1)","✓","✓", "✓", "") 
)

rows <- tibble(
  term = c("Dependent Variable", 
           #"Majority",
           #"President's Party",
           "All Legislators", 
           "Served At Least 2nd Term"),
  `(1)` = c("Count", "✓", ""),
  `(2)` =c("Count", "✓", ""),
  `(3)` = c("Count","", "✓"),
  `(4)` = c("Log(Count+1)", "✓", "") 
)

# subset
attr(rows, 'position') <- c(0, 20,21,22,23)

# call coefs
attr(rows, 'position') <- c(0, 24,25)


# HTML
modelsummary(models_total, stars = T, 
             add_rows = rows, 
             coef_map = cm,
             gof_omit = gof_omit,
             coef_omit = coef_omit,
             title = "The Effect of Expierence and Institutional Power on Legislator Contacts", 
             notes = list("Robust standard errors in parentheses, clustered by legislator.")) %>% row_spec(row = 1, bold = T, hline_after = TRUE) 
The Effect of Expierence and Institutional Power on Legislator Contacts
 (1)   (2)   (3)   (4)
Dependent Variable Count Count Count Log(Count+1)
Committee Chair 0.003 0.014 0.013 0.004
(0.006) (0.010) (0.010) (0.003)
Ranking Member −0.008 −0.006 −0.006 −0.002
(0.006) (0.009) (0.009) (0.003)
Prestige Committee −0.001 0.009 0.008 0.000
(0.003) (0.007) (0.007) (0.002)
First Year −0.001 −0.005 −0.004 0.000
(0.006) (0.011) (0.011) (0.005)
Second Year −0.006 −0.009 −0.008 −0.003
(0.006) (0.011) (0.011) (0.004)
Third Year 0.007 0.007 0.007 0.004
(0.007) (0.010) (0.010) (0.004)
Fourth Year −0.004 −0.004 −0.005 −0.002
(0.007) (0.010) (0.010) (0.003)
Fifth Year 0.014+ 0.014 0.015 0.003
(0.007) (0.009) (0.010) (0.003)
Sixth Year 0.000 −0.001 −0.001 −0.001
(0.007) (0.009) (0.009) (0.003)
Majority −0.003 −0.005 −0.005 −0.004*
(0.003) (0.005) (0.005) (0.002)
President’s Party −0.001 −0.001 −0.002 0.002
(0.003) (0.004) (0.004) (0.002)
All Legislators
Served At Least 2nd Term
Num.Obs. 412111 412111 388997 407272
FE: Year_x_Agency X X X X
FE: Legislator_x_Agency X X X
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001
Robust standard errors in parentheses, clustered by legislator.
# TeX 
modelsummary(models_total, # stars = T, 
             add_rows = rows,               
             coef_map = cm, 
             gof_omit = gof_omit, 
             coef_omit = coef_omit,
             output = "latex",
             #title = "\\\\footnotesize The Effect of Expierence and Institutional Power on Legislator Contacts", 
             notes = list("\\\\footnotesize Robust standard errors in parentheses, clustered by legislator.")) %>% 
  format_table() %>%
  write_lines(file = here::here("docs", "tables", "models_total.tex") )


beta <- m_total_dnd$coefficients %>%
  round(3) %>% 
  as_tibble(rownames = "beta") %>% 
  pivot_wider(names_from = beta)

se <- m_total_dnd$se %>%
  round(3) %>% as_tibble(rownames = "se") %>%
  pivot_wider(names_from = se)

Total Letters Predictions

  • Figures: figs/m-total-predicted-[1:4].png
# helper functions to plot predicted values
ggchair <- function(predicted = predicted) {
  
  predicted$chair %<>% str_replace("0", "Not Chair") %>% str_replace("1", "Chair")
  
  predicted$presidents_party %<>% str_replace("0", "Not President's Party") %>% str_replace("1", "President's Party")

predicted %>%  
  ggplot() + 
  aes(x = "", 
      y = estimate, 
      fill = chair,
      shape = chair,
      color = chair,
      ymin = conf.low,
      ymax = conf.high) + 
  geom_pointrange(position =  position_dodge(width = -.5)  )  + 
  scale_fill_viridis_d(begin = 0, end = .6, option = "cividis") +
  scale_color_viridis_d(begin = 0, end = .6, option = "cividis") +
      theme_minimal() + 
  theme(panel.border  = element_blank(),
        panel.grid.major.x = element_blank()) + facet_grid(. ~ presidents_party) 
}


ggtenure <- function(predicted = predicted) {
  
predicted %<>%
  # drop estimates from impossible values
      group_by(rowid, estimate, conf.low, conf.high, chair) %>% 
  mutate(sum = sum(first, second, third, fourth, fifth, sixth),
         more = ifelse(sum == 0, 1,0)) %>% 
    filter(sum < 2) %>% 
    pivot_longer(cols = c("first", "second", "third", "fourth", "fifth", "sixth", "more")) %>% 
    select(name, value) %>% 
    filter(value == 1) %>% 
    # clean up for presentation
    mutate(year_in_congress = name %>% 
             str_replace("more", "7\nor more") %>% 
             str_replace("sixth", "6") %>% 
             str_replace("fifth", "5") %>% 
             str_replace("fourth", "4") %>% 
             str_replace("third", "3") %>% 
             str_replace("second", "2") %>% 
             str_replace("first", "1") 
             ) %>% 
    ungroup()
  
# Ideally, we could plot against data, but there is so much variation that you can no longer distinguish differences in predicted values 
# predicted %<>% full_join(dcounts_tenure2 %>% mutate(predicted = NA))
  
predicted %<>% filter(chair == 0 | name %in% c("sixth", "more"))

  predicted$chair %<>% str_replace("0", "Not Chair") %>% str_replace("1", "Chair")
  

predicted %>%  
  ungroup() %>% 
  ggplot() + 
  aes(x = year_in_congress, 
      y = estimate, 
      shape = chair,
      color = chair,
      ymin = conf.low,
      ymax = conf.high) + 
  geom_pointrange(position =  position_dodge(width = -.3)  )  + 
  scale_fill_viridis_d(begin = 0, end = .6, option = "cividis") +
  scale_color_viridis_d(begin = 0, end = .6, option = "cividis") +
      theme_minimal() + 
  theme(panel.border  = element_blank(),
        panel.grid.major.x = element_blank()) 
  
}
predicted <- predictions(m_total_cross,
                         newdata = values)

# Cross-sectional predictions
predicted %>%
    mutate(estimate = estimate*90,
         conf.low = conf.low*90, 
           conf.high = conf.high*90) %>%
  ggchair() + 
  labs(title = "Predicted Total Letters per Year\n(Cross-Sectional)",
       x = "",
       y = "Predicted Total Letters Per Year",
       fill = "",
       color = "",
       shape = "") 

# by tenure
predicted <- predictions(m_total_cross,
                         newdata = values_tenure)

predicted %>%
    mutate(estimate = estimate*90,
         conf.low = conf.low*90, 
           conf.high = conf.high*90) %>% 
ggtenure() +
  labs(title = "Predicted Total Letters Per Year\n(Cross-Sectional)",
       x = "Years Serving in Congress",
       y = "Predicted Letters Per Year",
       fill = "",
       color = "",
       shape = "") 


# diff in diff 
predicted <- predictions(m_total_dnd,
                         newdata = values)


# predictions by chair and presidents' party
predicted %>%
    mutate(estimate = estimate*90,
           conf.low = conf.low*90, 
           conf.high = conf.high*90) %>%
  ggchair() + 
  labs(title = "Predicted Total Letters per Year\nDifference in Differences (Within-Legislator)",
       x = "",
       y = "Predicted Letters Per Year",
       fill = "",
       color = "",
       shape = "") 

# by tenure
predicted <- predictions(m_total_dnd,
                         newdata = values_tenure)

predicted %>%
    mutate(estimate = estimate*90,
           conf.low = conf.low*90, 
           conf.high = conf.high*90) %>%
ggtenure() +
    geom_line(aes(x = as.numeric(year_in_congress %>% str_sub(1,1)))) + 
  labs(title = "Predicted Total Letters Per Year\nDifference in Differences (Within-Legislator)",
       x = "Years Serving in Congress",
       y = "Predicted Letters Per Year",
       fill = "",
       color = "",
       shape = "") 
TotalTotalTotalTotal

Total

100 simulations

Subletting data to legislators in their first 6 years and setting the mean of the simulated DV (letters per legislator per year per agency) to their year (tenure) in office.

# A function to estimate diff in diff with a new random sample for DV 
sim_data <- dcounts_tenure %>% filter(tenure<7, tenure>0)

# sim_data$tenure

Cross sectional

simulation_cross <- function(x){
  
  set.seed(x)
  
  data <- sim_data %>% 
  # REPLACE DV WITH NORMAL DIST
  mutate(perYear = rnorm(n, tenure)) # year-2007))

m_total <-feglm (perYear ~ 
                       #chair + ranking_minority + prestige + 
                    first + second + third + fourth + fifth + sixth | Year_x_Agency, 
                    cluster = "Legislator",
                    # ALT vcov = hetero ~ ssc(cluster.adj = TRUE),
           data = data)

beta <- m_total$coefficients %>%
  #round(3) %>% 
  as_tibble(rownames = "beta") %>% 
  pivot_wider(names_from = beta)

return(beta)
}

# run 100 times 
sims <- map_dfr(.x = 1:100, .f =  simulation_cross)

# # plot estimates
# sims %>% 
#   ggplot() +
#   aes(x = first) + 
#   geom_density(fill = "grey")  #+  geom_vline(xintercept = 1)

sims %>% 
  ggplot() +
  aes(x = second) + 
  geom_density(fill = "grey") #+ geom_vline(xintercept = 2)

sims %>% 
  ggplot() +
  aes(x = third) + 
  geom_density(fill = "grey") #+ geom_vline(xintercept = 2)

sims %>% 
  ggplot() +
  aes(x = fourth) + 
  geom_density(fill = "grey") #+ geom_vline(xintercept = 2)

sims %>% 
  ggplot() +
  aes(x = fifth) + 
  geom_density(fill = "grey") #+ geom_vline(xintercept = 2)

sims %>% 
  ggplot() +
  aes(x = sixth) + 
  geom_density(fill = "grey") #+ geom_vline(xintercept = 2)

Diff in diff

simulation_dnd <- function(x){
  
  set.seed(x)
  
  data <- sim_data %>% 
  # REPLACE DV WITH NORMAL DIST
  mutate(perYear = rnorm(n, tenure)) # year-2007))

m_total_dnd <-feglm (perYear ~ 
                      # chair + ranking_minority + prestige + 
                      # majority + presidents_party 
                      first + second + third + fourth + fifth + sixth | Year_x_Agency + Legislator_x_Agency, 
                    cluster = "Legislator",
                    # ALT vcov = hetero ~ ssc(cluster.adj = TRUE),
           data = data)

beta <- m_total_dnd$coefficients %>%
  #round(3) %>% 
  as_tibble(rownames = "beta") %>% 
  pivot_wider(names_from = beta)

return(beta)
}

# run 100 times 
sims <- map_dfr(.x = 1:100, .f =  simulation_dnd)

# plot estimates
# sims %>% 
#   ggplot() +
#   aes(x = first) + 
#   geom_density(fill = "grey")  #+  geom_vline(xintercept = 1)

sims %>% 
  ggplot() +
  aes(x = second) + 
  geom_density(fill = "grey") #+ geom_vline(xintercept = 2)

sims %>% 
  ggplot() +
  aes(x = third) + 
  geom_density(fill = "grey") #+ geom_vline(xintercept = 2)

sims %>% 
  ggplot() +
  aes(x = fourth) + 
  geom_density(fill = "grey") #+ geom_vline(xintercept = 2)

sims %>% 
  ggplot() +
  aes(x = fifth) + 
  geom_density(fill = "grey") #+ geom_vline(xintercept = 2)

sims %>% 
  ggplot() +
  aes(x = sixth) + 
  geom_density(fill = "grey") #+ geom_vline(xintercept = 2)